The double well potential in quantum mechanics: a simple, numerically exact 

formulation 
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The double well potential is arguably one of the most important potentials in quantum mechanics, 
because the solution contains the notion of a state as a linear superposition of 'classical' states, a 
concept which has become very important in quantum information theory. It is therefore desirable 
to have solutions to simple double well potentials that are accessible to the undergraduate student. 
We describe a method for obtaining the numerically exact eigenenergies and eigenstates for such a 
model, along with the energies obtained through the Wentzel-Kramers-Brillouin (WKB) approxi- 
mation. The exact solution is accessible with elementary mathematics, though numerical solutions 
are required. We also find that the WKB approximation is remarkably accurate, not just for the 
ground state, but for the excited states as well. 



I. INTRODUCTION 

One of the most lasting impressions of quantum me- 
chanics in the public consciousness is that of the so-called 
Schrodinger cat paradox. 1 This notion, that physical sys- 
tems generally exist as a superposition of various states, 
lies at the heart of quantum mechanics, and is empha- 
sized, for example, when one studies the time evolu- 
tion of a system initially prepared in such a superpo- 
sition. Well publicized paradigms, such as the Einstein- 
Podolsky- Rosen (EPR) paradox, 2 , and the more modern 
tests of Bell's inequalities^^ are all manifestations of this 
superposition of states. 

A typical physics undergraduate student is likely to 
have been drawn to physics in the first place by pop- 
ular accounts of these and other phenomena. Part of 
maintaining this interest throughout their coursework is 
achieving the right balance between, on one hand, learn- 
ing the necessary mathematics and 'standard' solutions, 
as presented in most textbooks, and, on the other hand, 
actively engaging in problems with stimulating and vi- 
sual solutions that illustrate some of these exciting ideas. 
For example, Feynman used the ammonia molecule to il- 
lustrate the principle of superposition of states in quan- 
tum mechanics £ In systems like ammonia there are (typ- 
ically) two degenerate states; the system can reside in a 
superposition of both, and, most importantly, can tun- 
nel from one to the other. Tunnelling is omnipresent in 
quantum mechanics, and is the key ingredient in modern 
applications such as solid state devices (e.g. the diode), 
solar cells, and microscopes^ 

The example of the ammonia molecule is a rich prob- 
lem with many details one could focus on, but the pur- 
pose of this paper is to introduce the student to a more 
microscopic example. The simplest model for studying 
a superposition of nearly degenerate states is the one 
dimensional double well potential. This problem con- 
sists of a potential with two minima separated by a bar- 
rier. Model potentials of a double well are indeed treated 
in most textbooks, and even in many pedagogical-style 
articles 

The treatment of a double well potential in textbooks 
will vary from qualitative discussion (see, for example, 



Problem 2.47 in Griffiths^) to an analytical solution (for 
example, Section 8.5 in Merzbacher^ where, however, 
advanced mathematical steps are required) . More often, 
the Wentzel-Kramers-Brillouin (WKB) approximation is 
used (see, for example, Problem 8.15 in Griffiths^ and 
Problem 7.2, and 8.10 in Merzbacher^). In this paper 
we will present a method which allows students to nu- 
merically solve these tunneling problems for themselves 
without complicated mathematics. This method nec- 
essarily requires the students to access numerical proce- 
dures for diagonalizing a matrix, and has been explained 
in some detail in Ref. [1J]. The philosophy behind this 



teaching practice is that students can utilize familiar 
packages like Mathematica, Maple, and Matlab to both 
solve these problems and gain a visual and more mem- 
orable understanding of the solutions. This is currently 
not the case in most undergraduate classes, but asking 
students to carry out similar problems as outlined below 
will aid toward this goal. 

We begin by defining the problem in the next section, 
and outlining the numerically exact solution, following 
the treatment of the harmonic oscillator problem in Ref. 
[l4j |. The primary focus is the ground state splitting, 
but many other levels, including the wave functions, can 
be obtained. Our method of solution requires embed- 
ding the double well potential in an infinite well, whose 
presence has no effect on the ground state and low lying 
excited state solutions. We note that a variety of dou- 
ble well potentials can be studied this way; furthermore, 
there is no restriction on having an analytic form for the 
potential. This brings students to the point that very 
realistic potentials can be studied without further com- 
plication, bringing the classroom experience significantly 
closer to research problems. Indeed, one of the under- 
lying motivations for developing a set of problems like 
these, suitable to undergraduate training, is that these 
same methods are generally adopted in research practice. 

To make contact with more approximate analytic solu- 
tions, we outline the 'standard' WKB solution, defined 
by using connection formulas at linear turning points; 
here again, most authors have focused on the ground 
state splitting, but since we have numerically exact re- 
sults we will examine the WKB solution for all the states. 



This necessarily requires three separate calculations, one 
for the low-lying states inside the potential well, one 
for the higher states, where the energy of the state ex- 
ceeds the central barrier, and one for still higher energy 
states, where the vertical walls of the infinite well come 
into play. These latter states arise only because we have 
adopted the methodology of Ref. [3] , but the compar- 
ison is useful anyways for a proper assessment of the 
WKB approximation for all energy levels J£ This section 
is sufficiently detailed in part to illustrate that the WKB 
solution, when done fully, is actually very accurate. In 
the end, the WKB approximation, which is only valid 
for certain potentials, requires considerably more work 
than the numerically exact solution. 

The take-home message of this paper is that utiliz- 
ing the matrix mechanics of Ref. [14[ places the entire 
problem of the double well potential in full control of un- 
dergraduate students, and prepares them for even more 
complicated problems. 



II. THE DOUBLE WELL POTENTIAL: 
NUMERICALLY EXACT SOLUTION 

Several analytic forms for the double well potential are 
available for consideration, and are illustrated in Fig. 1. 
For a square double well the potential is defined by 

Vs , x) = {O if \x - (f - 6)| < | or \x - (f + 6)| < § 
| Vq otherwise, 

where the potential is centered around x — a/2, and the 
parameters 6 and c control the positioning and width of 
the individual wells, respectively. It is of course possi- 
ble to formulate this problem analytically, but it is te- 
dious, and more simply handled through the formalism 
described below. 

For curve-defined double wells we focus on two analyt- 
ical forms; the first one consists of two parabolic wells, 

^) = ^ E (k-f|-6) 2 , (2) 

where V cusp = ^muPb 2 is the value of the central max- 
imum, defined here in terms of an oscillator frequency, 
oj = y/k/m, where m is the mass of the particle, and k 
is the spring constant associated with the potential that 
exists on either side of the central maximum. The cen- 
tral maximum is located at x — a/2, and the minima 
are symmetrically located at x la ± = a/2 ±6, as is appar- 
ent in Fig. [T] The other potential, which arises in many 
symmetry breaking transitions, is a quartic function of 
position, 

n*0 = ^((*-f) 2 -6 2 ) 2 , (3) 

where V mayi is the value of the potential at the central 
maximum and a/2 ± 6 are the locations of the two local 



minima. In all cases the potential is symmetric about 
the point x = a/2, for reasons that will become clear 
below, and the two minima are centred around points 
a/2 ±6. 
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FIG. 1. (color online) The various possible double well po- 
tentials used in this problem are shown as a function of x. All 
are centred around x — a/2, and all have minima centered 
around x m — a/2 ± 6. The square double well potential is 
shown in red, with the analytic form given by Eq. fl]). The 
quadratic form is shown in blue, and is given by Eq. {2); 
it has a maximum (cusp) at the value V C us P . The quartic 
form is shown in green, and is described by Eq. ©; it has 
a smooth maximum at the value V^ax- All the potentials 
are embedded in an infinite square well potential, lying be- 
tween x = and x — a. This plot is schematic only, and 
the :r-axis is given in terms of the infinite square well width. 
In practice the potential would be provided with characteris- 
tic lengths, and then the infinite square well width would be 
chosen to ensure that states bound in the actual double well 
potential of interest would have negligible amplitude at the 
infinite wall boundaries (thus ensuring that the boundaries 
have no influence). 



In Ref. [14( we demonstrated that a number of po- 
tentials could be solved through matrix mechanics, by 
expressing the wave function in terms of a set of basis 
functions, and formulating an eigenvalue problem. Wc 
did this using a very simple and (perhaps) the most fa- 
miliar basis set known to undergraduate students at this 
stage, that of an infinite square well. Ref. [l4| solved the 
Schrodinger equation for familiar (but non-trivial) prob- 
lems, such as the harmonic oscillator potential. This was 
done so the reader could readily understand the effects 
of the confining walls of an infinite square well potential, 
along with the necessary truncation required for a matrix 
which is in principle infinite; however, these truncation 
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effects can be entirefy (and systematically) eliminated by 
experimentation. That is, for the low-lying states, the 
infinite square well width can be made sufficiently large 
so the solutions are no longer sensitive to the walls; fur- 
thermore, the dimension of the matrix can be increased 
until convergence to any desired accuracy is attained. 

In our formulation of the problem we require an in- 
finite square well; in Fig. 1 this has been chosen to lie 
between x = and x = a (also indicated with thick (red) 
vertical lines). The reason for placing this potential be- 
tween x — and x = a (rather than at —a/2 < x < a/2) 
is that the eigenstates of the infinite square well (without 
the double well potential) are more easily defined; they 
are conveniently given by 

Mx)={\H S ' m (^) [i0<X<a > (4) 
[ otherwise, 

with eigenvalues, 

The quantum number n — 1, 2, 3, . . . takes on a positive 
integer value, and we remind the reader that these states 
alternate between even and odd around the midpoint, 
x = a/2. 

Following Ref. [T3| , we expand the solution to our full 
problem in terms of this complete set: 

oo 

W) = ^2 Cm\lpm), (6) 
m=l 

and by inserting this into the Schrodinger equation and 
taking inner products with each basis state, we eventu- 
ally arrive at the eigenvalue equation 

oo 

^ ^ H nrn c m — Ec ni (7) 

m— 1 

where 

H nm =l^ n \H\ip m ) = 5 nm E^ 

+ - / dxsm )V(x)sm( (8) 

aj \ a I \ a I 

is the Hamiltonian matrix, and 5 nm is the usual Kro- 
necker delta function. 

The integration in Eq. © is elementary, regardless 
of which potential is used. Analytic expressions are pro- 
vided in the Appendix, but one should keep in mind that 
these integrals can also be readily obtained numerically. 
In fact, herein lies one of the advantages for students to 
learn this method of solving the Schrodinger equation, as 
they can in principle solve for many realistic potentials. 

Note that, just like in the harmonic oscillator problem, 
H nrn is symmetric and peaked along the matrix diago- 
nal. One might worry about the effects of truncation 
for this matrix problem; however, for large n, the diag- 
onal elements increase as n 2 . The off-diagonal elements 



quickly become negligible in comparison, and for large 
n the energy levels approach those of an infinite square 
well, with an offset which accounts for the presence of 
the potential (the offset is equal to the average of the 
potential being considered). 

Equation ([7]), with the use of equation (|A3|) . (|A6D . 
or (|A12[) in the Appendix can be evaluated up to some 
cutoff for given parameters of the potential (for exam- 
ple, the parameters hui/E^ and b/a in the case of the 
quadratic potential - Eq. ©). Equations (IA3p . (|A6I) . 
and (IA12|) (each pertaining to a particular potential) 
form the elements of the matrix which is fed into an 
eigenvalue/eigenvector solver (see for example, Numeri- 
cal Recipes^ or defined functions in software packages 
such as Matlab, Mathematica, or Maple). Typical re- 
sults are shown in Fig. [2] for each of the potentials. To 
arrive at such results, students will have to truncate the 
matrix at some cutoff, N; a cutoff as low as N — 20 leads 
to inaccuracies of less than 0.2% for the ground state en- 
ergies of the potentials ([U, and (O, but the square well 
potential (JXJ) requires a larger cutoff of N = 50 to be 
just as accurate. The reason for this is that the double 
square well potential presents a much narrower structure 
compared to the other two; therefore an accurate solu- 
tion requires contributions from basis states with higher 
frequency components. These higher frequency states 
are only available at larger quantum number, n. By in- 
creasing the matrix size we can systematically produce 
the exact results to arbitrary precision. 

There is a lot of room for experimentation; if one 
elects to study the two square wells, for example, stu- 
dents will discover that the convergence rate decreases 
dramatically as the well widths become narrower. This 
occurs for the same reason as it did in the case of a sin- 
gle square well^ namely that higher energy basis states 
are necessary not because they have higher energy, but, 
as explained above, because they oscillate in real space 
with shorter and shorter periods, and these more refined 
spatial variations are required to describe wave functions 
that are confined to narrower wells. 

As one increases the maximum matrix size, N, the 
matrix solution will provide more eigenvalues and eigen- 
vectors. It is important for students to realize that large 
eigenvalues will not correspond exclusively to a double 
well potential, but will also include the effects of the 
infinite square well in which the former is embedded. 
This means the higher eigenvalues and eigenvectors will 
closely resemble those pertaining to the infinite square 
well, as students can readily verify. Corrections to the 
infinite square well results correspond to accounting for 
the average of the double well potential, a result that 
students can readily verify for each potential. 

Returning to Fig. [21 we have chosen the potential pa- 
rameters to obtain essentially the same ground state 
eigenvalue for each potential (see caption for details) . In 
all cases the tunneling barrier that results is sufficiently 
large that the first excited state will be almost degener- 
ate (slightly higher) with the ground state energy. We 
show the ground state and first excited state wave func- 
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FIG. 2. (color online) The ground state (solid black curves) and first excited state (purple circles) wavefunctions for the 
three potentials (shown as dashed red, blue, and green curves, for the square, quadratic, and quartic double well potentials, 
respectively) are shown. All potentials are symmetric about x = a/2 with minima located at x m = a/2 ± 6, just as in Fig. [T] 
A value of b/a — 0.2 is used so that the infinite walls embedding the double well have no effect on the lower eigenstates; 
meanwhile, the height of the tunnelling barrier is such that several excited states occur in doubles (see Fig. 3 below). The 
values of Vb, V CU s P , and U ma x were determined by requiring the same ground state energy for all three wells. The square barrier 
potential (red curve in Fig. 2a) has an extra degree of freedom, the individual well widths, c. We used c/a — 0.1 to make these 
results similar to those for the other two potentials. With a ground state energy, E\/E^ — 50, b/a = 0.2, (and c/a = 0.1 for 
the square well), we determined V /e[ 0) « 253.6, V cusp /e[ 0) sa 987.0, and V m ^/E[ 0) « 260.2. Note that the wave functions 
for the ground state and the first excited state are essentially identical, except for the differing symmetry. Also, the three sets 
of wave functions are very similar to one another. The slightly narrower square double well result (a) is explained in the text. 



tions in each case. It is clear that they are generically 
the same; the ground state has even parity, and shows 
a large amplitude to be in either the left or right well. 
The first excited state has the same characteristics, ex- 
cept that it has odd parity. These results illustrate that 
the characteristics for low-lying states are qualitatively 
similar, independent of the precise form of the double 
well potential. Upon closer inspection it is also clear 
that the individual peaks of the wave function for the 
square double well potential (Fig. 2a) are slightly nar- 
rower than the other two. This is consistent with the fact 
that this well is more confining than the other two, and 
also with the fact noted earlier, that more basis states 
were required to obtain this solution accurately. 

As mentioned above, one can continue with any (or all 
three) potentials and examine more characteristics of the 
tunneling between the two sides, or of the excited states. 
In the next section we focus on one particular potential, 
the one which is quadratic in x, given by Eq. and ex- 
amine the eigenvalues more closely. Results for this po- 
tential are shown in Fig. [3] Focusing on just the numeri- 
cal results, note that as one moves up in energy, we pass 
through the different regimes; at low energies the levels 



are split and come in almost degenerate pairs, while at 
higher energies the levels spread out, corresponding to 
a regime where the central region is no longer forbid- 
den. Finally, at very high energies the spectrum begins 
to approach quadratic (in n) behaviour, corresponding 
to an infinite square well potential. In the next section 
we further discuss these results after providing a brief 
summary of the WKB results. 



III. THE DOUBLE WELL POTENTIAL: WKB 
APPROXIMATION 

As indicated above, there are three distinct regions in 
this problem; these regions are not completely obvious 
from the numerical solution, but their presence needs to 
be explicitly accounted for when undertaking a WKB 
calculation. The reader should note that the WKB so- 
lution for a double well potential is presented in most 
textsjii but usually the treatment is confined to a cal- 
culation of the splitting between the two lowest energy 
levels. This (tiny) splitting is connected to the tunneling 
probability between the left and right side wells, giving 
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rise to issues analogous to those that arise in the case of 
the Schrodinger cat J£ The WKB approximation for all 
three regions differs in details - depending on whether 
the energy level is below the central maximum (region 
1), E < V(a/2) = V^usp, above the central maximum 
but below the value of the double well potential at the 
walls (region II), V cusp < E < V(a), or above this value 
(region III), E > V(a). 

It is worthwhile, given that numerically exact solutions 
are so readily available with this method, to investigate 
how accurate the WKB solution really is. The solutions 
physically relevant to a double well all lie in region I; 
beyond this our method of solution (i.e. the existence of 
infinite walls) affects the numerical solutions^ Nonethe- 
less once one works out the WKB solution in region I, 
the superb agreement with the numerically exact results 
encouraged us to investigate the degree of agreement in 
the other two regions. As we see below, the agreement 
is excellent there as well, and we therefore describe each 
of these in the subsections below. These calculations 
are slightly more involved; the undergraduate student 
may want to follow the solution in region I, while a full 
survey of all three regions may be more appropriate for 
students working on specialized projects (related to the 
WKB approximation), or for graduate students. 



A. region I 

Region I (below the central maximum) is the region 
treated in most textbooks and pedagogical papers, since 
the splitting in energy levels that results from tunnelling 
is the most important outcome of the double well poten- 
tial. While textbooks often deal with the ground state, 17 
we wish to describe not just the ground state, but all 
eigenstates whose energy lies in Region I. Thus, following 
the methodology of Ref. [TtJ , but avoiding their approx- 
imations, leads to the following transcendental equation, 



E n '± = huj( 



where Q = Q/E^ for all energies, and the quantity 
is itself a function of the energy that we are seeking: 
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n'± 



(10) 



Eqs. (|9ll0p must be solved iteratively. Note that there 
is an integer index n' which is associated with two 
(split) levels, and we reserve the label n (which appears 
in Fig. 3, for example) to enumerate all the eigenval- 
ues. In practice the iteration converges very quickly, 
as in the ground state, for example, E « V(a/2), so 

(j) ~ 2V f^ 2 ^ ! as derived in Ref. [I?} , and Eq. ([9]) provides 
a closed form solution for the two energies. Similarly, the 



argument of the inverse tangent function is very small, so 
we obtain for the ground state, in the case of a relatively 
high central barrier, 



A^E + -E^^e- 



where 



2V(a/2) 



7T 

T 



E 



(0) 



(11) 



(12) 



all of which is well known^ 

The numerical results (iterated solutions to Eqs. 
(|9I10I) ) are shown in Fig. 3 as green cross-hairs (inside 
the squares in Region I), and clearly coincide with the 
(numerically) exact results (the centres of the squares). 
Of course the splitting is a much finer detail and will 
be examined further below. Note that the ID quan- 
tum harmonic oscillator result is also shown (dashed pink 
line starting near the origin); the numerically exact and 
WKB results nicely follow that of a single harmonic os- 
cillator for most of Region I. 



B. region II 

Region II involves intermediate energies, as indicated 
in Fig. 3. The WKB procedure is straightforward, and 
results in the transcendental equation: 



E n — 



fe»(n-|) 



l + § 



V(§) 



1 - 



E„ 



sm 



E„ 



(13) 



where the 'bar' over each energy indicates it is normal- 
ized to the infinite square well ground state energy, e[ ^ . 
Once again, this equation requires iteration, but we find 
convergence is readily attained within a handful of it- 
erations. The results are also plotted in Fig. 3 as pink 
asterisks (in the squares in Region II). In this region 
there is a distinct departure from the result describing a 
single harmonic oscillator, but the WKB results are in 
very good agreement with the numerically exact ones. 



C. region III 

Region III consists of energy levels above the point 
where the double well potential is cut off by the walls. 
This region exists only because of our method of nu- 
merically exact solution, but we include a comparison 
nonetheless. The WKB approximation is now specific to 
that involving two vertical walls. The equation govern- 
ing this result is 



E n = 



[l + *i(a) + « 1 (f) + &(a) + &(§)]' 



(14) 



5 





I 


I 


i i i 


4000 


□ Exact 


hw = 
b/a = 


100 
0.2 


3000 






x WKB III 


F 


region III 






c n 








2000 


—I 

region II 


* WKB II j 


i — lijd * 


1000 


1 







region I „ *lsgp 


+ WKB I 







J>l?l 

I 


i 


i i i 



10 20 30 40 50 60 

n 



FIG. 3. (color online) Eigenvalues obtained by the numerical diagonalization of the double well potential ([2]) embedded in 
an infinite square well (see Fig. [TJ. The numerically exact results are shown with red squares (the actual points are at the 
center of the squares). We used a truncation size of JV = 200 to obtain these results, but a much smaller size would suffice as 
well, especially for the low lying eigenvalues. We also show with various symbols the WKB approximations appropriate to the 
three different regions, as discussed in the text (they are also color coded here). They are remarkably accurate on this energy 
scale. Also shown is a dashed (pink) straight line starting from near the origin; it is given by E n — (hQ/2)(n — 1/2), which 
represents the result expected for a single harmonic oscillator potential. This expectation is relevant when the two parts of 
the double well are 'very far away' from one another, and so the two parts of the double well appear to behave independently. 
It is clear that the levels in Region I satisfy this criterion very well. (Note: the factor of 2 (in to/2) is required because there 
are two wells, and for the same reason we need twice as many quantum numbers. One should think of the pairs of values in 
region I as sharing a single quantum number which is the average of the two.) Also shown as the black double-dotted curve 
is the analytic expression, given by Eq. (| ltjp ; this was derived through first order perturbation theory and is expected to be 
very accurate for the higher levels in Region III. Indeed, it is already fairly accurate at n — 60. 



where the corrections labelled Si, i — 1,2, tend to be 
small, but are dependent on the energy being sought. 
These corrections are given by 




Eqs. ([Tlf and (fT5|) are iterated to convergence with very 
little difficulty. These results are shown as the (blue) 
i's (in the squares in Region III), and are in excellent 



agreement with the numerically exact results. In fact, 
a good starting approximation is obtained by using first 
order perturbation theory (with the double well potential 
as the perturbation and particle in the infinite square 
well as the unperturbed Hamiltonian), and one obtains 14 

This curve is also shown in Fig. |3] as a double-dotted 
(black) curve, and, though not as accurate as the WKB 
approximation, works fairly well for levels in Region III. 
These results indicate the remarkable accuracy of the 



G 



WKB approximation, not just for the ground state, but 
for the excited states as well. The latter is not too sur- 
prising, since the WKB is a semiclassical approximation, 
and ought to improve as the quantum number increases; 
however, this is not often emphasized in comparisons of 
this sort. 



D. The energy splitting 

While the comparison of the WKB approximation to 
the numerically exact results in Fig. 3 may look impres- 
sive, the scale there is fairly coarse; indeed a 'large' er- 
ror in the splitting of the levels in Region I would go 
unnoticed on such a plot. Therefore we take a closer 
look at this energy splitting, since, ultimately, it is this 
splitting which plays a primary role in the tunnelling 
phenomenon. 

The WKB approximation motivates a suitable defini- 
tion for the energy splitting observed for the low lying 
states. Using Eq. f) as a guide, we define the energy 
splitting to be A„< = E n > + — where the index 

n' refers to the pair of eigenvalues clearly discernible in 
Fig. [3j A comparison of these energy splittings with the 
WKB approximation spans several decades, so we find 
it more convenient to instead focus on the exponent <j) 
that appears in Eq. ©. Within the WKB approxima- 
tion this is defined by the iterated result of Eq. (|10[) (in 
combination with Eq. ©), and we will henceforth refer 
to this as </vwkb- On the other hand, the exact result 
is defined by 
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FIG. 4. (color online) Plot of the exponent <f> using the WKB 
approximation (Eqs. © and (|10p - symbols) and as defined 
by the numerically exact results (Eq. (|17p - colored curves). 
The WKB approximation is in excellent agreement with the 
numerically exact solutions for the four pairs of split states, 
even down to the point where the barrier is about to dis- 
appear for the 4th pair of states (lower left). Also shown 
are analytical approximations for the WKB approximation as 
discussed in the text; these are plotted as thin solid (black) 
curves, and are given by Eq. (|18[) . The double-dotted (black) 
curve includes a further correction discussed below Eq. (|18[l , 
and is included only for the 4th pair of states, as it is not 
needed for the others. These analytic approximations to the 
WKB iterated solution work very well for all the pairs of 
states shown. 



(17) 



where the E n i± are obtained numerically. The numer- 
ically exact results are compared with the WKB in 
Fig. |31 the exponent cf> is shown plotted as a function 
of 2V(a/2)/(7icj) = (tt 2 /2)hu){b/ a) 2 . This is the barrier 
height with respect to the approximate ground state en- 
ergy level. In Fig. 0] this ratio is changed by varying 
lu such that 40 < Tluj < 100. Clearly the agreement is 
(very) excellent over the entire range of hui shown, for 
the first 4 pairs of energy levels. Note that this is true 
even at the lower range of x (i.e. uj just below the central 
maximum), where the 4th pair is just about to enter re- 
gion II, so this is the limit of essentially no barrier at all. 
Closer scrutiny shows that the WKB approximation very 
slightly overestimates the level of 4> and therefore slightly 
underestimates the energy splitting, compared to the nu- 
merically exact result. Moreover, we obtain a very good 
estimate for </vwkb by expanding Eq. (TTU)) assuming 
that E n ,± k, hui(n' - 1/2), and that x = 2V{a/2)/{Tuui) 
(the barrier height in Fig. 4) is large; the result is 







H) 






\2n- 1 J 



and these are also shown in Fig. 4 (solid black curves). 
This expression works very well for the ground state, 
but begins to deteriorate for the excited states, espe- 
cially at the lower range of x (w near the central maxi- 
mum) . Expanding to the next order adds a term *' 2 " 8 ~ 1 ' ) 
to Eq. ([T8|) . which is shown for the 4th pair only as a 
black double-dotted curve — this correction essentially 
restores the approximate value to the fully iterated WKB 
value and even more so for the other three pairs (not 
shown). We conclude that the WKB approximation is 
indeed very accurate, even when examining finer energy 
scales associated with the splitting of the low-lying lev- 
els that results from a very small tunneling amplitude 
between the two wells. 



IV. SUMMARY 

In this paper we have provided an illustration of the 
numerical solution to a non-trivial problem, the double 
well potential. The focus has not been on the numerical 
implementation; on the contrary, we have used this ex- 



7 



ample in the classroom, with the expectation that stu- 
dents have software tools (Mathematica, Maple, Mat- 
Lab, etc.) readily available. In practice this is not al- 
ways the case, but our belief is that this expectation is 
slowly becoming an integral part of an undergraduate 
education in science or engineering. Assuming this is 
so, then this method of solution requires a conceptual 
understanding of mathematical constructs that are typ- 
ically taught in the first year of university: integrals in 
calculus and linear algebra. We emphasize that concep- 
tual understanding only is required, since implementa- 
tion of both of these concepts is readily accomplished 
through software tools. The alternative is to first un- 
derstand the (conceptual) solution to differential equa- 
tions, and then utilize numerical software to implement 
these solutions; this is a much more difficult task, since 
an understanding of differential equations lags behind 
introductory calculus and linear algebra by about two 
years. In our experience with teaching^ Quantum Me- 
chanics, much of the students' understanding of differ- 
ential equations emerges from advanced physics courses 
(such as QM). This in turn focuses attention away from 
the physical systems and onto understanding the math- 
ematics, which is precisely what we hope to avoid with 
the present approach. 

We have emphasized the characteristics that arise in 
the solution of the double well potential that one can 
teach conceptually through perturbation or variational 
procedures. However, by using the (numerically exact) 
method presented here, the student gains several advan- 
tages. First, the solution is numerically exact, so linger- 
ing doubt about the validity of the approximate solutions 
is no longer possible. Second, this method of solution can 
be applied to any problem of this kind, so the poten- 
tial doesn't even have to be analytically defined. Third, 
this method of solution is an implementation of what is 
usually taught formally in most undergraduate quantum 
mechanics courses, and certainly reinforces the more ab- 
stract concepts of Hilbert space, etc. Finally, this is pre- 
cisely the method of solution used in research (at least 
in condensed matter problems, one rarely solves a dif- 
ferential equation — usually problems are formulated in 
terms of matrices) . Therefore practice with these meth- 
ods can better equip a student for research in physics. 

It is clear from the double well examples presented, 
that for the low-lying states, the solutions have much in 
common. The wave function amplitude is peaked around 
each of the wells, so classically speaking, the particle can 
be viewed as being in a superposition of two states; in 
the first the particle is in the left well, while in the sec- 
ond the particle is in the right well. With respect to 
these solutions, it is clear by visual inspection of Fig. [21 
for example, that the vertical walls that constitute the 
embedding infinite square well potential have no influ- 
ence on the solutions, as all wave functions shown have 
decayed to zero at the walls. 

In addition we have carried out a detailed analysis 
of the WKB approximation for a quadratic double well 
potential, as this method is often used for such prob- 



lems in real-world applications. We have determined the 
WKB energy levels self-consistently, and found that this 
approximation is remarkably accurate, not just for the 
ground state, but for the excited states as well. 
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Appendix A: matrix elements 

As remarked in the text, since it is assumed that the 
student is using a software package to diagonalize the 
matrix, it is likely straightforward to utilize the same 
package to evaluate the matrix elements either analyti- 
cally or numerically. Alternatively, to proceed 'by hand,' 
we use the identity 



2 sin a sin b = cos (a — b) — cos (a + b) , 



(Al) 



to simplify the integral in equation ([5]). 

For the square double well, the potential is non-zero 
and constant (Vq) over the intervals 



<x< - - b 



a , c a 

2- b+ 2 <X< 2 

a c 

- + b+ - <x< a, 



b- 



(A2) 



so it is a simply a matter of integrating cos (x) over these 
intervals. The result is (with E± as the unit of energy) 



(ft. 



{Hs) n jE? ] = 5 n 



V ( 2c 2(-l)" (1irnb\ (imc\ 
' ' 1 1 sin I I 



(0) 



a irn 



sin 

a I \ a J 



+ (1 - S nm ) ^ -^y (q(n + m) - q(n - m)j , (A3) 



where 



1 / im\ ( nnb\ ( irnc 

q(n) = — cos — cos sin 

y ' n \ 2 J \ a I 2a 



(A4) 



We have used the subscript 's' to denote the 'square' 
double well potential. 

For the quadratic and quartic potentials only integrals 
with even powers of x times cos x are required: 



,1/2 

Ij (n) = dx x J cos (nirx) , 

Jo 



(A5) 
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where the integration need only extend to the half-way 
point because of the symmetry of the potential. The 
result for the quadratic potential is 



(h 2 



_ (H 2 ) T 



r(0) 



= 8„ 



+ (1 - 8 nm 
where 
fn{x) = 1 

and 



E\ 

1 + (-1)"+™ 



48 V E[ 



Tiuj 



AO) 



f I--- 

In{ 2 a 



(A6) 



6(-l) 
(nn)'- 



6x 



,(-!)' 



(7m) 2 



12x 2 , 



(A7) 



1 /(_l)("-™)/2 (•_ 1 ^(n+m)/2 

2 \ (n — my [n + my 



(__1)(«-™)/ 2 _ 1 (_ 



+m)/2 



(n — m) 2 



(n + m) 2 



(A8) 



We can formulate an equivalent expression to equa- 
tion (|A6I) , where if n and m are both odd, 



{hi 



{H2) n J E i 



(0) 



V CUSP 8mn[l-g(2+(^ + ^)(-l)^)] 



(0) 



[<57t(to — n)(m + n)] ' 



(A9) 



where <5 = b/a, and if n and m are both even, 
(M„ m = {H2) n jE { ? ] = 

V cusp 8mrc[l-2<5(l-(-l)^)] 



(0) 



E 



and when n = 



[Sn(m — n)(m + n)] ' 



(A10) 



1 / 1 



(^)nn/^i 
(-1)' 



(0) 



cusp 



(nTT) 



4 0) 
S 2 \ 12 



2(n7r) 2 



(All) 



and for all other n and to, (h 2 ) nm = 0. 

Here we have used the subscript '2' to denote the 
quadratic potential. For equation (| A6|) we have adopted 



the (dimensionless) energy ratio, -757, to denote the 

E i 

strength of the potential (recall V cusp = \muj 2 b 2 , where 
V cusp is the value of the central maximum and the loca- 
tions of the minima at x m — a/2 ± b defines the parame- 
ter b). Also note that, just like in the harmonic oscillator 
problem, the g nrn (|A8|) remain of order unity close to the 
diagonal. Nonetheless, as in that problem, for large n, 
the diagonal elements grow as n 2 , so the off-diagonal el- 
ements become negligible in comparison, and the energy 
levels approach those of an infinite square well, with an 



offset to account for the double well potential. This is 
done by averaging the potential and adding the aver- 
age to the upper eigenvalues of the infinite square well, 
thus producing shifted infinite well energies for the upper 
eigenvalues of the embedded double well potential. 

Finally, for the quartic potential ([3]), using the iden- 
tity (jAll) and the substitution y = x/a — 1/2 results in 
the expression, 

(M„ m - m nJ E[ 0) = S nm n 2 + cos {-±+1!^ 

x li % [(-!) m ^ - "0 - J(n + to)] , (A12) 

where 5 = b/a, and 

J(n) = h{n) - 25 2 I 2 {n) + (5 4 / (n), (A13) 

where Ij(n) is defined in equation (|A5|) . These integrals 
are readily obtained through integration by parts, inte- 
gral tables, or analytic integral solvers such as Maple 
and Mathematica. The required expressions are 

[ I / 7f ji \ 

Io(n) = -S n . + (1 - <5„. )— sin — , (A14) 
2 nn V 2 / 



h{n) 



24 

/_L_ 2 



Snfi + (1 - <5«,o) X 



\47rn (7m) 3 



sin 



/7m 

v~2~ 



(7m) 2 



, (A15) 



and 



/4 ( n ) = WW + ( x ~ S n,o) x 
loO 

1 3 



24 



I67™ (7m) 3 (rrn) 5 
1 12 \ /mi 



2(7rn) 2 (im) 4 



cos 



V 2 



(A16) 



Alternatively, the following equivalent expression can 
be obtained using analytic integration solvers such Maple 
and Mathematica (with appropriate simplifying assump- 
tions such as n and to being integers): 



(0) 



1 fi 



1 



1/1 



S 2 \6 (mr) 2 ) 5 4 \80 4(mr) 2 2(nn) 



1 



(0) 



1 



+ (1 - S nm ) 



v m 



E\ 



(o) 



2mn (l + (-l)" +m ) 
\tt5 (n — to) (n + to)] 2 



7T 2 (n 2 - to 2 ) 2 - 48(n 2 + to 2 ) 
\tt6 (n — to) (n + to)] 2 



(A17) 



As stated in the text, these can be readily checked by 
a simple numerical evaluation of the original integral. 
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